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Abstract 

A thermostat of the Nose-Hoover type, based on relative velocities and a local definition of 
the temperature, is presented. The thermostat is momentum-conserving and Galilean- invariant, 
which should make it suitable for use in Dissipative Particle Dynamics simulations, as well as 
nonequilibrium molecular dynamics simulations. 



* Invited Talk at Symposium on Progress and Future Prospects in Molecular Dynamics Simulation - In 
Memory of Professor Sliuiclii Nose 



I. INTRODUCTION 

The original papers of Nose provided a new perspective on the generation of sta- 

tistical ensembles by dynamical simulation. They showed that a deterministic set of equa- 
tions of motion, involving just one or two extra degrees of freedom, can sample configura- 
iions from the canonical ensemble. This complements the stochastic method of Andersen 
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which generates the canonical ensemble by periodic reselection of velocities from the 
Maxwell-Boltzmann distribution. The work of Hoover P, further clarified the nature of 
the isothermal dynamical equations and how to derive them (see also the contribution of 
Hoover to these proceedings). The Nose-Hoover equations incorporate a dynamical friction 
coefficient, whose fluctuations are driven by the difference between the instantaneous ki- 
netic temperature (defined through the sum of squares of particle velocities) and the desired 
temperature. 

This paper presents a new thermostat of the Nose-Hoover type, based on an instanta- 
neous temperature which is calculated as a weighted sum of squares of relative velocities of 
atom pairs. The frictional term in the equations of motion also enters in pairwise fashion, 
conserving momentum, and making the dynamics invariant to a Galilean transformation of 
velocities. There are two areas in which such a thermostat may be useful. The first area 
is nonequilibrium molecular dynamics (NEMD). With a conventional Nose- Hoover thermo- 
stat, it is necessary to apply the friction term to the peculiar velocities of the particles, i.e. 
the difference between the true velocities and the local streaming velocity of the fluid at 
the particle positions. Failure to do this can lead to unphysical (thermostat-induced) be- 
haviour |g] such as the stabilisation of string phases 0, Q or the generation of steady-state 
antisymmetric stress . One solution to these problems is to use a proflle- unbiased thermo- 
stat, which requires a self-consistent determination of the streaming velocity in the course 
of the simulation 
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. A thermostat based on the relative velocities of nearby pairs of 
atoms may avoid, or at least ameliorate, the problem. The second area of application of 
the pairwise thermostat is dissipative particle dynamics (DPD). Here, in order to preserve 
hydrodynamic behaviour, it is essential for any thermostat to conserve momentum, and a 
pairwise form is one way of achieving this. This paper concentrates on the DPD case, since 
the suggestion of usi ng a pairwise Nose-Hoover thermostat was flrst made in this context by 
Stoyanov and Groot |lO]. However, it should be borne in mind that the thermostat may be 
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applied equally well to, e.g. Lennard- Jones fluids. 

The paper is organised as follows. Section |n] contains a brief summary of DPD, concen- 
trating on the temperature control aspects. Section lllll derives the equations of motion for 
the pairwise Nose-Hoover thermostat, and also summarizes the equations for a thermostat 



based on the conflgurational temperature, due to Braga and Travis for comparison. 
Section IIVI presents the results of some preliminary tests for DPD simulations. Finally, 
section IVl contains the conclusions. 



II. DISSIPATIVE PARTICLE DYNAMICS 



DPD [12, |13| has become a popular tool for simulating the behaviour of both simple and 
complex fluids. It consists of the solution of the classical equations of motion for a system 
of interacting particles, together with a set of stochastic and dissipative forces which control 
the temperature and allow one to choose the viscosity. For a simple fluid the equations may 
be written H Q, Q 

ri = Vi= Pi/rrii (la) 
P^ = Mr) - e^(r,p) + aR,{r,p) , (lb) 

where r and p stand for the complete set of coordinates and momenta. The so-called conser- 
vative forces fi are derived from a pair-potential term in the Hamiltonian = —{dH/dvi) 
and so may be written as fi = fij^ with fji = —fij. In DPD these pair forces usually 
take the form 

fij = awij = aw{rij) , with w{r) = w{r)r (2a) 

{( 1 — r /rc) r < Tc 
(2b) 
r > Tc 

Here rij = Vi — Vj, r = \r\, r = r/r. The parameter a determines the strength of the 
conservative interactions, and rc is the cutoff. 

The dissipative forces — ^ Vi are also written in pairwise fashion Vi = X^j^i with 
Vji = —Vij, usually defined thus: 

V^j = {vij ■ Wij)wij = w{rijf{v,j ■ rij)rij (3) 



where Vij = Vi — Vj. A choice has been made here to use the same weighting function w{r) 
as in the specification of conservative forces. aRi is short for the random "forces", which 
also act between pairs, with a weight function w{r); the strength parameter a is related 
through the fiuctuation-dissipation theorem to the friction coefficient ^ and the temperature 
k^T (see for more details). The pairwise nature of all these forces guarantees the 

momentum conservation necessary to ensure hydrodynamic behaviour: in other words, the 
dynamics is Galilean- invariant. The particles represent fiuid regions, rather than individual 
atoms and molecules: the softness and simplicity of the interactions permit the use of a long 
time step, compared with conventional molecular dynamics. This, and the acceleration of 
physical processes compared with those seen in more realistic simulations, gives an advantage 
of several orders of magnitude, at the cost of a very rough mapping onto specific molecular 
properties. 

A slightly more general view of DPD treats it as conventional molecular dynamics using 
soft potentials, supplemented by a momentum-conserving thermostat which acts between 
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pairs. Lowe [15] takes this approach, rather than solving the above equations. Instead, each 
timestep At involves the following operations. 

1. Positions and momenta are advanced using = pi/rrii, Pi = fi. 

2. Every pair ij (in random order, and possibly subject to a distance dependent weight 
or range function) is examined and, with probability P = uAt, the momenta are 
updated: Pi : = pi + Apij, pj ■ = Pj — Apij, with 



Apij = rrii 



C\ kBT/rriij - {vij ■ ri 



where ( is picked from a Gaussian distribution with zero mean and unit variance, and 
rriij = niimj/ {rrii + rrij). 

This procedure periodically reselects the component of the relative velocity along from the 
Maxwell-Boltzmann distribution corresponding to reduced mass rriij. The key parameter is 
the stochastic randomization frequency u: high values of u give effective temperature control, 
but also a high viscosity; low values give very weak temperature control while allowing the 
viscosity to^be low. The thermostat is closely related to the one originally proposed by 
Andersen 



4 



Recently, Stoyanov and Groot have proposed a modification of tlie above metliod: 
tlie fraction (1 — P) of pairs wfiich do not have their relative velocities stochastically updated, 
are instead thermalized by a deterministic method. For each such pair, a dissipative force is 
calculated and used to correct the momenta during the deterministic part of the step, incor- 
porating a temperature-dependent controlling factor. Finally, the Lowe velocity reselection 
process is applied to the remaining fraction P of pairs as before. The idea of Stoyanov and 
Groot is to give more control over the separate effects of thermalization, namely temperature 
control and changing viscosity. Stoyanov and Groot jl^ call the deterministic part of their 
thermostat "Nose- Hoover" , but actually it is not of this form, and has not been shown to 
generate the canonical ensemble. It may be noted that an algorithm resembling that of Nose 
and Hoover was also described by Besold and Mouritsen |16[. 

III. PAIRWISE NOSE-HOOVER THERMOSTAT 
A. Derivation of Equations of Motion 

The purpose of this paper is to present a Galilean-invariant thermostat of the Nose- 
Hoover type, which generates the canonical ensemble. The derivation is a straightforward 
implementation of the approach of Hoover P|, anc a special case of the generalized Nose- 



Hoover equations discussed by Kusnezov et al. 
assumed to be of the form 



and Martyna et al. ^^J. The result is 



ri = pi/mi (4a) 
P^ = Mr)-mr,p) (4b) 
i = G^ir,p) (4c) 

with the Vi{r,p) given by eqn 0. Eqns ()4a|) and ()4b|) are written down by analogy with 
eqns (^. The random forces are dropped, the friction coefficient ^ is now an additional 
dynamical variable, and the right-hand side of eqn (jlc|) is the object of the derivation. 
This is obtained from the generalized Liouville equation for the (stationary) phase space 
distribution function g{r,p,^) 

i:|--(^^.)+i:|--(^p.)4(p?)-o. (5) 
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The ansatz is made that G^{r,p) in eqn (Pc|) depends only on positions and momenta, so 
d^/d^ = 0. Direct substitution shows that equation (0) is satisfied by the product form 

p(r,p,0 oc exp{-H{r,p)/kBT}exp{-lQ^e/kBT} 

where is an arbitrary constant, provided 

i ^ 

= 5Z • - ikBT/m,j)w{njf) 

i j<i 

= 5Z 5Z ^('"^i)^ [{'"ii ■ ^iiT ~ kBT/niij] . (6) 

i j<i 

Once more, the reduced mass rriij appears. The term in square brackets vanishes if an aver- 
age is taken over the canonical momentum distribution. The equation has a straightforward 
physical interpretation, acting to damp the difference between the instantaneous tempera- 
ture corresponding to the component of relative velocity Vij along the inter-particle vector, 
and the canonical ensemble average of this quantity. The prefactor controls the "thermal 
inertia" in the same way as the corresponding parameter in the conventional Nose-Hoover 
method, and the function w gives a higher weighting to closer pairs. There is a conserved 
"energy function" 

Hti{r,p,i,^i,) = H{r,p) + ig^^^ + if^ where if^ = ik^iT^wirijf /mij (7) 
as may be checked by time differentiation and direct substitution of the equations of motion. 



B. Integration Algorithm 

It is not the aim here to discuss the optimal algorithm for integration of the equations of 
motion (Q. Instead, the simplest modified velocity-Verlet algorithm that is commonly 
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used in DPD |16[, is adopted: 



Pi■.= pr.= P^ + ^At{fi-^V,) P^(^At) (8a) 

e:=^=e + iAtG^ e(^At) (8b) 

Vi := ri + Atpi/mi ri{At) (8c) 

/. := Mr) MAt) (8d) 

Vr.= V,{r,p) ViiAt) (8e) 

G^:=G^{r,p) G^{At) (8f) 

p, := + iAt (/, - ei^) p,(At) (8g) 

^:=^+lAtG^ e(At) (8h) 



Steps (lHe|l - (j81ij) may be iterated to convergence, because the momenta at time t + At should 
be used in the evaluation of and Vi. However, because of the expense of calculating the 
pairwise terms, in DPD it is usual to stop after one evaluation of the expressions above, and 
this is the approach adopted here. Some might prefer a strictly reversible integrator |2ol |. 
while others favour the Runge-Kutta method: consideration of these possibilities is deferred. 



C. Configurational Nose-Hoover Thermostat 



The canonical ensemble result 



dU 



(9) 



has been known for many years 



temperature in simulation 



211 and has recently been used to define a configurational 
22 , 23 1 and experiment 24 , l^l • Recently, one of us Q| has 



suggested monitoring this quantity as an indicator of lack of equilibrium due to excessive 
timesteps in DPD. It is natural to consider applying a thermostat to control this variable 
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27l | and here the equations of motion of Braga and Travis ll| are used: 



Vi = Pi/rrii + fifi{r) 
Pi = Mr) 
A = G,{r) 



(10a) 
(10b) 
(10c) 
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where 



G, 




d dU 



dvj dvj 



(11) 



The quantity /i plays the role of a fluctuating mobility: that is, a proportionality between 
force andjirift velocity, as seen in the "position Langevin equation" or Schmoluchowski 



equation 



ll|. Once again there is a conserved "energy function" 



H^{r, p, /i, (f^) = H{r, p) + \Qf,^? + where 0^ = ^ik^T ^ 



d dU 



(12) 



Braga and Travis (1]| have presented a simple integration algorithm for these equations, 
which we use here. The canonical distribution may also be shown to be a steady-state 
solution of the above equations of motion, and they share with the thermostatted equations 
of section UlIAl the property of Galilean invariance. 



IV. RESULTS 

Tests have been carried out using the standard "water" DPD model Q]: the potential 
strength parameter in eqn Q was set to a = 25, with simulation units deflned so that 
m = 1, k-^T = 1, Tc = 1. A system of = 250 particles was simulated in cubic periodic 
boundaries. Timesteps in the range 0.005 < At < 0.06 were used, with run lengths up to 
1000 reduced time units. For the pairwise Nose-Hoover thermostat, inertia parameters in 
the range 0.2 < Q^/N < 8.0 were studied. For the conflgurational Nose-Hoover thermostat, 
inertia parameters in the range 2000 < Q^/N < 80000 were used. 

These thermostats allow one to check the accuracy of the integration by monitoring the 
conserved energy-function eqns ((Tj) and (fT^. Figure [T] shows that, at timesteps At > 0.02 
(very conservative by DPD standards), there is a signiflcant drift in this quantity: the rate 
of increase is roughly proportional to At^ at large At. This problem has been noted before 
by Hafskjold et al. |23], and it is not associated with the thermostatting, because the same 
behaviour is seen using the simple velocity Verlet algorithm. The cause seems to be the 
relatively strong discontinuity in force derivatives at the cutoff of the DPD potential iboj. In 
DPD, and in MD with a thermostat, this tends to be camouflaged. The present thermostats 
perform as well as (in fact, slightly better than) velocity Verlet in this respect. 

The oscillation of the internal energy of the particles (potential plus kinetic) reflects the 
flow of energy into and out of the thermal reservoir, and this is influenced by the choice of 

8 



FIG. 1: Rate of change of energy-like function as a function of timestep At, plotted on log-log 
scales. Circles: pairwise Nose- Hoover thermostat with inertia parameter Q^/N = 0.8. Squares: 
configurational Nose-Hoover thermostat with inertia parameter Q^/N = 2 x 10^. Diamonds: ve- 
locity Verlet algorithm, with no thermostatting. The lines correspond to Af^ power law behaviour. 




At 

thermal inertia parameter. Typical results are shown in Fig. |21 and they show the expected 
behaviour. There are damped oscillations: the runs lengths employed here are typically long 
compared with the relaxation rate, while the timesteps are small enough to cope with the 
oscillations. In this range, the precise choice of thermal inertia is not critical. 

The simulation-averaged values of kinetic temperature (defined through the total ki- 
netic energy) and configurational temperature (defined by eqn Q) are shown in Fig. 01 
When the kinetic temperature is controlled, lack of equilibrium is indicated by the configu- 
rational temperature, which increases by as much as 10% at the largest timesteps studied. 
These results simply confirm what has been seen before j^] : a measured kinetic temperature 
close to the desired value should not be taken as a guarantee that the system is at equilib- 
rium. The form of the increase in may be understood semi-quantitatively by considering 
the effect of non-zero-timestep velocity- Verlet dynamics on the phase portrait of a simple 
harmonic oscillator ^26^. Conversely, when the configurational thermostat is imposed, the 
measured kinetic temperature is significantly reduced when the timestep is too large. This 
effect may be understood in a similar way by considering harmonic oscillator velocity- Verlet 
dynamics: for a given positional amplitude, the momentum amplitude is reduced as the 
timestep increases. When both thermostats are applied together, not surprisingly, both T^. 
and Tk are controlled well, up to the highest timesteps studied. This deserves further inves- 
tigation, but it would be over-optimistic to suppose that the other degrees of freedom in the 
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FIG. 2: Oscillation of internal energy E = H as a function of time t. Upper panel: pairwise 
Nose-Hoover thermostat with inertia parameter: Q^/N = 0.2 (solid line); Q^/N = 0.8 (dashed 
line); Q^/N = 2.0 (dot-dashed line. Lower panel: configurational Nose-Hoover thermostat with 
inertia parameter: Qf,/N = 8 x 10^ (solid hne); Q^/N = 2 x 10^ (dashed line); Q^/N = 4 x 10^ 
(dot-dashed line). 
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ication to complex fluids, simulations of the same lipid bilayer model 
3fl[ 0| have been carried out with the new thermostat. Here, the 



system are at equilibrium. 

To illustrate the app 
studied previously j2f 

solvent water is represented as before, and each lipid molecule has the form of a 7-bead chain 
HTg in which a-repulsion parameters between hydrophilic "head" beads (H), hydrophobic 
"tail" beads (T), and "water" beads (W) are chosen to produce the desired behaviour js^. 
Harmonic bond-stretching potentials, and angle-bending potentials, act within the lipid 
molecules. The measured temperatures of the different types of DPD bead are shown in 
Figure m The results are consistent with those obtained before [2^ and show how dangerous 
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FIG. 3: Kinetic temperature Tk (open symbols) and configurational temperature Tc (filled symbols) 
as functions of timestep At for three different therniostatting regimes: (a) pairwise Nose-Hoover 
thermostat with inertia parameter Q^/N = 0.4; (b) configurational Nose-Hoover thermostat with 
inertia parameter Q^/N = 4 x 10^; (c) both thermostats simultaneously. 
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FIG. 4: Kinetic temperature (open symbols) and configurational temperature (filled symbols) 
as functions of timestep At for membrane simulations using three different thermostatting regimes: 
(a) pairwise Nose-Hoover thermostat with inertia parameter Q^/N = 0.4; (b) configurational Nose- 
Hoover thermostat with inertia parameter Q^/N = 4 x 10^; (c) both thermostats simultaneously. 
Different symbols represent different DPD particle types: circles, H, Tg; squares, Ti, T5; diamonds, 
T2, T3, T4; triangles, water. 
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it is to rely on thermostats to equilibrate the system when the timestep is too large: the 
different bead types have significantly different kinetic and configurational temperatures in 
all cases. Actually, for this simple model, the cause of the problem, and the remedy, are 
well understood. The intramolecular potentials within the lipid chains are too strong to be 
handled by the longer timesteps; this problem is easily addressed by using multiple timestep 
methods j32]. However, this example serves to illustrate possible pitfalls which may occur 
in the general case. 



V. DISCUSSION 

The derivation of Section IIIII establishes the canonical ensemble as a stationary distri- 
bution for the coordinates and momenta subject to the equations of motion although 
it does not prove that it is unique, nor guarantee that a system will converge towards this 
distribution 12 , [3| • The result is easily generalized to apply to a subset of pair inter- 
action, dimply by =ettmg = fa the omitted pa.., m^ng thi= suitable to combine 
with the Lowe method as envisaged by Stoyanov and Groot [10[ . (Interestingly, Ref . [5| con- 
tains exercises on incorporating a weighting factor, and on considering a subset of degrees 
of freedom, for the conventional Nose- Hoover thermostat). Nose-Hoover chains may easily 

n 

be added to further control the dynamics Il8j |. 

The preliminary results presented above indicate that the pairwise Nose-Hoover thermo- 
stat behaves as should be expected, and may be useful in both DPD and conventional MD 
/ NEMD simulations. A feature of the proposed thermostat, shared by the configurational- 
temperature thermostat, is the absence of peculiar velocities: this may provide a more 
satisfactory way of controlling the temperature than the conventional Nose-Hoover ther- 
mostat in the case of fluid flows, since only local relative velocities are used to define an 
instantaneous temperature. 
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